I Предварительный анализ данных

if (interactive() && !str_ends(getwd(), "R/Statistical_Analysis_2022/CITY_US/Tkachenko")) {
    setwd("R/Statistical_Analysis_2022/CITY_US/Tkachenko")
}

data <- read_excel("../CITY_US.STD/CITY_shortname.xls")
data[data == "NA"] <- NA
data[, -(1:2)] <- data.frame(lapply(data[, -(1:2)], as.numeric))
fullnames <- names(read_excel("../CITY_US.STD/CITY.xls"))

1 Разобраться в том, что означают признаки.

print(fullnames)
##  [1] "CITY City"                                                                            
##  [2] "STATE State"                                                                          
##  [3] "AREA Land area, 1990 (sq.miles)"                                                      
##  [4] "R_AREA Land area, 1990 (sq.miles), ranks"                                             
##  [5] "POP92 Population,1992"                                                                
##  [6] "R_POP92 Population,1992, ranks"                                                       
##  [7] "POP80 Population,1980"                                                                
##  [8] "R_POP80 Population,1980, ranks"                                                       
##  [9] "POP_CH Population growth rate,1980-1992"                                              
## [10] "R_POP_CH Population growth rate,1980-1992, ranks"                                     
## [11] "POPDEN Population per square mile, 1992"                                              
## [12] "R_POPDEN Population per square mile, 1992, ranks"                                     
## [13] "OLD % older 65 years"                                                                 
## [14] "R_OLD % older 65 years, ranks"                                                        
## [15] "BLACK Black population,1990"                                                          
## [16] "R_BLACK Black population,1990, ranks"                                                 
## [17] "BLACK% % Black population,1990"                                                       
## [18] "R_BLACK% % Black population,1990, ranks"                                              
## [19] "ASIAN Asian or Pacific Islander population,1990"                                      
## [20] "R_ASIAN Asian or Pacific Islander population,1990, ranks"                             
## [21] "ASIAN% % Asian or Pacific Islander population,1990"                                   
## [22] "R_ASIAN% % Asian or Pacific Islander population,1990, ranks"                          
## [23] "HISP Hispanic population, 1990"                                                       
## [24] "R_HISP Hispanic population, 1990, ranks"                                              
## [25] "HISP% % Hispanic population, 1990"                                                    
## [26] "R_HISP% % Hispanic population, 1990, ranks"                                           
## [27] "BORN_F % persons of foreign born"                                                     
## [28] "R_BORN_F % persons of foreign born, ranks"                                            
## [29] "LANG % of persons, speaking language other than English at home, 1990"                
## [30] "R_LANG % of persons, speaking language other than English at home, 1990, ranks"       
## [31] "HH1 % 1-person households, 1990"                                                      
## [32] "R_HH1 % 1-person households, 1990, ranks"                                             
## [33] "FAMIL1 % one-parent family households, 1990"                                          
## [34] "R_FAMIL1 % one-parent family households, 1990, ranks"                                 
## [35] "DEATH Infant death rate per 1000 live births,1988"                                    
## [36] "R_DEATH Infant death rate per 1000 live births,1988, ranks"                           
## [37] "CRIME Crime rate per 100000 population, 1991"                                         
## [38] "R_CRIME Crime rate per 100000 population, 1991, rate"                                 
## [39] "SCHOOL % of elementary and high school enrollment in public schools, 1990"            
## [40] "R_SCHOOL % of elementary and high school enrollment in public schools, 1990, ranks"   
## [41] "DEGREE % of adults with a bachelor's degree or higher, 1990"                          
## [42] "R_DEGREE % of adults with a bachelor's degree or higher, 1990, ranks"                 
## [43] "INCOME Median household income, 1989"                                                 
## [44] "R_INCOME Median household income, 1989, ranks"                                        
## [45] "ASSIST % of households, receiving public assistance,1989"                             
## [46] "R_ASSIST % of households, receiving public assistance,1989, ranks"                    
## [47] "POVERT % of persons below poverty level, 1989"                                        
## [48] "R_POVERT % of persons below poverty level, 1989, ranks"                               
## [49] "OLB_BIL % of housing units built 1939 or earlier, 1990"                               
## [50] "R_OLD_BI % of housing units built 1939 or earlier, 1990, ranks"                       
## [51] "OWNER Median value of specified owner-occupied housing units ($), 1990"               
## [52] "R_OWNER Median value of specified owner-occupied housing units ($), 1990, ranks"      
## [53] "RENTER % renter-occupied housing units,1990"                                          
## [54] "R_RENTER % renter-occupied housing units,1990, ranks"                                 
## [55] "GROSS Median gross rent of specified renter-occupied housing units ($), 1990"         
## [56] "R_GROSS Median gross rent of specified renter-occupied housing units ($), 1990, ranks"
## [57] "CONDOM % condominimum occupied housing units, 1990"                                   
## [58] "R_CONDOM % condominimum occupied housing units, 1990, ranks"                          
## [59] "TRANSP % of workers using public transportation"                                      
## [60] "R_TRANSP % of workers using public transportation, ranks"                             
## [61] "UNEMP Unemployment rate, 1991"                                                        
## [62] "R_UNEMP Unemployment rate, 1991, ranks"                                               
## [63] "LABOR Labor force - % change 1980-1990"                                               
## [64] "R_LABOR Labor force - % change 1980-1990, ranks"                                      
## [65] "LAB_F Female civilian labor force participation rate, 1990"                           
## [66] "R_LAB_F Female civilian labor force participation rate, 1990, ranks"                  
## [67] "MANLAB % of civilian labor force emploied in manufacturing, 1990"                     
## [68] "R_MANLAB % of civilian labor force emploied in manufacturing, 1990, ranks"            
## [69] "TAXE City government taxes per capita, 1991"                                          
## [70] "R_TAXE City government taxes per capita, 1991, ranks"                                 
## [71] "TEMPER Average daily temperature in July"                                             
## [72] "R_TEMPER Average daily temperature in July, ranks"                                    
## [73] "PRECEP Annual precipitation"                                                          
## [74] "R_PRECEP Annual precipitation"

2 Отобрать признаки

names_interesting <- c("AREA", "POP80", "POP92", "POPDEN", "CRIME", "BORN_F", "POVERT", "INCOME", "UNEMP", "TEMPER")

data <- data %>% select(all_of(c("CITY", "STATE", names_interesting)))

print(head(data))
## # A tibble: 6 × 12
##   CITY  STATE  AREA  POP80  POP92 POPDEN CRIME BORN_F POVERT INCOME UNEMP TEMPER
##   <chr> <chr> <dbl>  <dbl>  <dbl>  <dbl> <dbl>  <dbl>  <dbl>  <dbl> <dbl>  <dbl>
## 1 NEW_… NY     309. 7.07e6 7.31e6  23671  9236   28.4   19.3  29823   8.6   76.8
## 2 LOS_… CA     469. 2.97e6 3.49e6   7436  9730   38.4   18.9  30925   9     74.3
## 3 CHIC… IL     227. 3.01e6 2.77e6  12185    NA   16.9   21.6  26301   8.4   75.1
## 4 HOUS… TX     540. 1.60e6 1.69e6   3131 10824   17.8   20.7  26261   6.1   83.5
## 5 PHIL… PA     135. 1.69e6 1.55e6  11492  6835    6.6   20.3  24603   8     76.7
## 6 SAN_… CA     324  8.76e5 1.15e6   3546  8537   20.9   13.4     10   6.2   71

3 Определить вид признаков

Город и штат качественные, остальные количественные, ранги были порядковыми.

find_mode_freq <- function(x) {
    x <- x[!is.na(x)]
    return(max(tabulate(match(x, x))))
}

print(data %>% summarise(across(
    all_of(names_interesting),
    find_mode_freq
)))
## # A tibble: 1 × 10
##    AREA POP80 POP92 POPDEN CRIME BORN_F POVERT INCOME UNEMP TEMPER
##   <int> <int> <int>  <int> <int>  <int>  <int>  <int> <int>  <int>
## 1     1     1     1      2     1      4      2      1     5      3
print(sort(data$UNEMP))
##  [1]  2.3  3.2  3.8  3.9  4.1  4.4  4.6  4.7  4.7  4.8  4.8  5.0  5.0  5.0  5.0
## [16]  5.1  5.1  5.1  5.3  5.4  5.4  5.4  5.4  5.4  5.6  5.7  5.9  5.9  5.9  6.0
## [31]  6.0  6.1  6.1  6.1  6.1  6.2  6.2  6.4  6.4  6.5  6.7  6.7  6.7  6.8  6.9
## [46]  6.9  7.0  7.1  7.2  7.2  7.3  7.3  7.3  7.5  7.6  7.7  7.7  7.7  8.0  8.1
## [61]  8.4  8.4  8.5  8.6  9.0  9.0  9.4  9.5  9.6 10.0 10.4 10.6 10.7 11.0 11.9
## [76] 12.6 13.1

Все количественные буду считать непрерывными. возможно UNEMP непрерывный с плохой точностью.

4 не актуально

5 Построить matrix plot

if (interactive()) pdf("ggpairs_unedited.pdf")
ggpairs(
    data[, -(1:2)],
    lower = list(continuous = wrap("points", alpha = 0.5, size = 0.3)),
    diag = list(continuous = "barDiag")
)

if (interactive()) dev.off()

7 outliers

Убираю outliers:

Помечаю некорректные данные в INCOME как NA. Удаляю город из Аляски за плотность населения. Флорида выделсется на BORN_F-INCOME Гаваи выделяются низким уровнем безработицы. странно?

data$INCOME[data$INCOME < 100] <- NA
data <- data %>% filter(STATE != "AK")

6 Несимметричные распределения

Функция, которая логарифмирует, если это сделает выборку симметричнее

log_asymmetric <- function(x) {
    if (skewness(x, na.rm = TRUE) < abs(skewness(log(x), na.rm = TRUE))) {
        print("default")
        return(x)
    } else {
        print("logged")
        return(log(x))
    }
}

Автоматически логарифмирую то что имеет асимметрию и длинный хвост справа

data_logged <- data %>%
    mutate(across(all_of(names_interesting), log_asymmetric))
## [1] "logged"
## [1] "logged"
## [1] "logged"
## [1] "logged"
## [1] "logged"
## [1] "logged"
## [1] "default"
## [1] "default"
## [1] "logged"
## [1] "default"
if (interactive()) pdf("ggpairs_logged.pdf")
ggpairs(
    data_logged[, -(1:2)],
    lower = list(continuous = wrap("points", alpha = 0.5, size = 0.3)),
    diag = list(continuous = "barDiag")
)

if (interactive()) dev.off()

8 однородность

выглядит однородно.

9 не актуально

10 всякие характеристики

print_characteristics <- function(x) {
    list(
        mean = mean(x, na.rm = TRUE),
        var = var(x, na.rm = TRUE),
        skewness = skewness(x, na.rm = TRUE),
        kurtosis = kurtosis(x, na.rm = TRUE)
    )
}

sapply(data_logged %>% select(all_of
(names_interesting)), print_characteristics)
##          AREA      POP80     POP92     POPDEN    CRIME      BORN_F    POVERT   
## mean     4.754638  12.9069   13.01226  8.257586  9.206558   1.959747  18.11842 
## var      0.6658984 0.5142076 0.4464288 0.5099096 0.06798296 0.8435754 34.68872 
## skewness 0.1317646 1.453879  1.743773  0.1015321 0.1476689  0.3081197 0.2243883
## kurtosis 2.609264  6.033026  6.885752  2.846858  3.088721   2.256848  2.646043 
##          INCOME     UNEMP      TEMPER    
## mean     25282.48   1.873674   77.60132  
## var      14321837   0.09967048 37.59853  
## skewness -0.2196584 -0.2186434 -0.1426746
## kurtosis 2.353211   3.659464   3.747206

II О виде распределений и о сравнении распределений

1 выполнить первое второе задание в логичном порядке

2 анализ вида распределения признаков

2. Описать распределения признаков визуально (с помощью нормальной бумаги и PP-plot) и по критерию хи-квадрат и другим критериям.

Сюда входит: normal probability plot (что это такое?), проверка по критериям Лиллиефорса, AD, хи-квадрат, Шапиро-Уилка. По критерию хи-квадрат, а также визуально по PP-plot можно проверить и гипотезы о согласии с другими распределениями, например, логнормальным.

Рассматриваем нормальность

все плотности.

if (interactive()) pdf(file = "all_density_plots.pdf")
ggplot(melt(data_logged, id = c("CITY", "STATE")), aes(value)) +
    theme_bw() +
    geom_histogram(aes(y = ..density..), bins = 15) +
    # geom_density() +
    facet_wrap(~variable, scales = "free") +
    geom_line(aes(y = dnorm(value,
        mean = tapply(value, variable, mean, na.rm = TRUE)[PANEL],
        sd = tapply(value, variable, sd, na.rm = TRUE)[PANEL]
    )), color = "red") +
    theme(legend.position = "bottom") +
    labs(x = "", y = "")

if (interactive()) dev.off()

Normal probability plot

if (interactive()) pdf(file = "normal_probability_plot.pdf")
ggplot(melt(data_logged, id = c("CITY", "STATE")), aes(sample = value)) +
    stat_qq_point(size = 2) +
    geom_abline() +
    facet_wrap(~variable, scales = "free")

if (interactive()) dev.off()

PP-plot

if (interactive()) pdf(file = "PP_plot.pdf")
ggplot(melt(data_logged, id = c("CITY", "STATE")), aes(sample = value)) +
    stat_pp_point(size = 2) +
    facet_wrap(~variable, scales = "free") +
    stat_pp_line()

if (interactive()) dev.off()

Всякие тесты на нормальность

test_all <- function(data, test, ...) {
    print(test(1:10, ...)$method)
    print(data %>% summarise(across(
        everything(),
        function(x) {
            test(x, ...)$p.value
        }
    )))
}

test_all(data_logged[, -(1:2)], shapiro.test)
## [1] "Shapiro-Wilk normality test"
## # A tibble: 1 × 10
##    AREA     POP80       POP92 POPDEN CRIME BORN_F POVERT INCOME UNEMP TEMPER
##   <dbl>     <dbl>       <dbl>  <dbl> <dbl>  <dbl>  <dbl>  <dbl> <dbl>  <dbl>
## 1 0.757 0.0000113 0.000000310  0.673 0.725 0.0892  0.675  0.463 0.487  0.368
test_all(data_logged[, -(1:2)], lillie.test)
## [1] "Lilliefors (Kolmogorov-Smirnov) normality test"
## # A tibble: 1 × 10
##    AREA  POP80    POP92 POPDEN CRIME BORN_F POVERT INCOME UNEMP TEMPER
##   <dbl>  <dbl>    <dbl>  <dbl> <dbl>  <dbl>  <dbl>  <dbl> <dbl>  <dbl>
## 1 0.536 0.0118 0.000633  0.291 0.580  0.148  0.646  0.866 0.746  0.766
test_all(data_logged[, -(1:2)], ad.test)
## [1] "Anderson-Darling normality test"
## # A tibble: 1 × 10
##    AREA    POP80       POP92 POPDEN CRIME BORN_F POVERT INCOME UNEMP TEMPER
##   <dbl>    <dbl>       <dbl>  <dbl> <dbl>  <dbl>  <dbl>  <dbl> <dbl>  <dbl>
## 1 0.492 0.000167 0.000000752  0.370 0.620 0.0362  0.830  0.692 0.697  0.371
test_all(data_logged[, -(1:2)], pearson.test)
## [1] "Pearson chi-square normality test"
## # A tibble: 1 × 10
##    AREA  POP80  POP92 POPDEN CRIME BORN_F POVERT INCOME UNEMP TEMPER
##   <dbl>  <dbl>  <dbl>  <dbl> <dbl>  <dbl>  <dbl>  <dbl> <dbl>  <dbl>
## 1 0.697 0.0377 0.0179  0.664 0.728  0.265  0.729  0.857 0.897 0.0626

итак, все кроме населения распределено нормально (или логнормально)

1. Описать разницу между численностью населения в 1980 и 1992 годах визуально (с помощью ящиков с усами) и по критериям. Добавить признак, который разбивает города на “холодные” и “жаркие”. Сравнить их по характеристикам.

добавим признак

temper_split <- mean(data_logged$TEMPER)
data_logged$F_TEMPER <- factor(ifelse(data_logged$TEMPER < temper_split, "cold", "hot"))

3 Сначала имеет смысл посмотреть на сравнение распределений в группах с помощью ящиков с усами

if (interactive()) pdf("boxplot_pop.pdf")
ggplot(
    melt(data_logged, id.vars = "F_TEMPER", measure.vars = c("POP80", "POP92")),
    aes(x = variable, y = value, color = F_TEMPER)
) +
    geom_boxplot() +
    scale_color_manual(values = c("blue", "red"))

if (interactive()) dev.off()


if (interactive()) pdf("boxplot_pop_all.pdf")
ggplot(
    melt(data_logged, id.vars = "F_TEMPER", measure.vars = c("AREA", "POP80", "POPDEN", "CRIME", "BORN_F", "POVERT", "INCOME", "UNEMP", "TEMPER")),
    aes(x = variable, y = value, color = F_TEMPER)
) +
    geom_boxplot() +
    facet_wrap(~variable, scales = "free") +
    scale_color_manual(values = c("blue", "red"))
## Warning: Removed 11 rows containing non-finite values (`stat_boxplot()`).

if (interactive()) dev.off()

Население растет, правый ящик с усами выше.

4 начинать нужно с t-критерия…

leveneTest(value ~ variable, data = melt(data_logged[c("POP80", "POP92")]), center = "mean")
## No id variables; using all as measure variables
## Levene's Test for Homogeneity of Variance (center = "mean")
##        Df F value Pr(>F)
## group   1  0.3172 0.5741
##       150
t.test(
    data_logged$POP80,
    data_logged$POP92,
    var.equal = TRUE,
    paired = TRUE
)
## 
##  Paired t-test
## 
## data:  data_logged$POP80 and data_logged$POP92
## t = -4.7135, df = 75, p-value = 1.098e-05
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -0.14988791 -0.06083053
## sample estimates:
## mean difference 
##      -0.1053592

Мы не доверяем критерию фишера потому что он работает только для нормальных. применим тест левена, он не отвергает гипотезу что дисперсии равны, как и ожидалось для почти одинаковых распредлений. считаем что дисперсии равны и применяем t-test.

5 если бы выборка было маленькая или дисперсии были неравны, пришлось бы применять непараметрический t-test Вилксона.

wilcox.test(data_logged$POP80, data_logged$POP92, paired = TRUE)
## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  data_logged$POP80 and data_logged$POP92
## V = 730, p-value = 0.0001492
## alternative hypothesis: true location shift is not equal to 0

6 Смотрите на результаты применения критерия Манна-Уитни

Манна-Уитни = Вилксона, этому тесту не нужны дополонительные требования, но он менее мощен чем t-test.

7 критерий Колмогорова-Смирнова

еще есть тест колмогорова смирнова, который сравнивает распеределения.

ks.test(data_logged$POP80, data_logged$POP92)
## 
##  Exact two-sample Kolmogorov-Smirnov test
## 
## data:  data_logged$POP80 and data_logged$POP92
## D = 0.15789, p-value = 0.3012
## alternative hypothesis: two-sided

он не отвергает.

8 сравнение зависимых выборок

это все и было сравнение зависимых выборок

из-за того что тесты парные, они смогли отвергнуть гипотезы о том, что средние равны, без этого гипотезы не отвеграются. статистика зависимого критерия больше статистики независимого (при положительной корреляции), поэтому мощнее.

3 Об анализе зависимостей

1 вспомним pairs plot

if (interactive()) pdf("ggpairs_logged.pdf")
ggpairs(
    data_logged[, -(1:2)],
    lower = list(continuous = wrap("points", alpha = 0.5, size = 0.3)),
    diag = list(continuous = "barDiag")
)

if (interactive()) dev.off()
if (interactive()) pdf("ggpairs_logged_colored.pdf")
ggpairs(
    data_logged[, -(1:2)],
    aes(color = F_TEMPER),
    lower = list(continuous = wrap("points", alpha = 0.5, size = 0.6)),
    diag = list(continuous = "barDiag"),
    columns = names_interesting
) +
    scale_color_manual(values = c("blue", "red")) +
    scale_fill_manual(values = c("blue", "red"))

if (interactive()) dev.off()

Матрица корреляций пирсона

data_logged <- data_logged %>%
    relocate(TEMPER, AREA, POP92, POP80, POPDEN, BORN_F, UNEMP, POVERT, CRIME, INCOME)


if (interactive()) pdf(file = "cor_matrix_pearson.pdf")
ggplot(melt(cor(data_logged %>% select(-CITY, -STATE, -F_TEMPER), method = "pearson", use = "pairwise.complete.obs"))) +
    geom_raster(aes(x = Var2, y = Var1, fill = value)) +
    geom_text(aes(x = Var2, y = Var1, label = round(value, 2))) +
    scale_fill_gradient2() +
    theme_dark() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

if (interactive()) dev.off()

Матрица корреляций спирмана

if (interactive()) pdf(file = "cor_matrix_spearman.pdf")
ggplot(melt(cor(data_logged %>% select(-CITY, -STATE, -F_TEMPER), method = "spearman", use = "pairwise.complete.obs"))) +
    geom_raster(aes(x = Var2, y = Var1, fill = value)) +
    geom_text(aes(x = Var2, y = Var1, label = round(value, 2))) +
    scale_fill_gradient2() +
    theme_dark() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

if (interactive()) dev.off()
if (interactive()) pdf(file = "all_density_groupped_plots.pdf")
ggplot(
    melt(data_logged, id = c("CITY", "STATE", "F_TEMPER")),
    aes(value, color = F_TEMPER)
) +
    theme_bw() +
    # geom_histogram(bins = 15) +
    geom_density() +
    facet_wrap(~variable, scales = "free") +
    geom_line(aes(y = dnorm(value,
        mean = tapply(value, variable, mean, na.rm = TRUE)[PANEL],
        sd = tapply(value, variable, sd, na.rm = TRUE)[PANEL]
    )), color = "black", linetype = 2) +
    theme(legend.position = "bottom") +
    labs(x = "", y = "") +
    scale_color_manual(values = c("blue", "red"))

if (interactive()) dev.off()

data_logged$F_TEMPER <- NULL

# plot_all_density(data_logged_melt, aes(value, color = F_TEMPER)) +
#     scale_color_manual(values = c("blue", "red"))
cor(data_logged$POP80, data_logged$POP92, method = "pearson")
## [1] 0.9628711
cor(data_logged$POP80, data_logged$POP92, method = "spearman")
## [1] 0.9155161